{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Python Practice 11"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This Python Practice will be dedicated to non-linear curve fitting using the `lmfit` library. After installing it, we will do a short introduction to the library and then we will fit a simple function using it."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Prerequisites"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's first install the `lmfit` library. You can do this directly from Jupyter Notebook by running the following command:\n",
    "\n",
    "```\n",
    "!pip install lmfit\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip install lmfit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import lmfit as lm\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Theoretical background"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `lmfit` package provides simple tools to help you build complex fitting models for non-linear least-squares problems and apply these models to real data. The way it works, is that for any data you want to fit, you need to define a model function, which has some free parameters to be optimized. The `lmfit` library will then use the Levenberg-Marquardt algorithm to minimize the residuals between the data and the model. The total fitting error can be expressed in a chi-square statistic:\n",
    "\n",
    "$$ \\chi^2 = \\sum_{i=1}^{N} \\left( \\frac{y_i - f(x_i)}{\\sigma_i} \\right)^2 $$\n",
    "\n",
    "where $y_i$ is the data, $f(x_i)$ is the model function, and $\\sigma_i$ is the uncertainty in the data. The goal is to minimize this statistic by adjusting the free parameters of the model function."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Therefore, while fitting the data, you need to follow these steps:\n",
    "\n",
    "1. Define the model function with free parameters.\n",
    "2. Estimate the uncertainties in the data.\n",
    "3. Define the residual function.\n",
    "4. Define the initial guesses for the free parameters.\n",
    "5. Minimalize the residuals w.r.t. the free parameters."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Defining the data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this Python Practice we will generate synthetic data of a simple sinusoidal function with some noise. The function we will use is:\n",
    "\n",
    "$$ y = 3\\sin(5 x + \\pi/3)$$\n",
    "\n",
    "We only assume that we know in advance that the function is a sinusoidal one, but we don't know the amplitude, frequency, and phase. Therefore, our model function will have three free parameters: `a`, `b`, and `c`:\n",
    "\n",
    "$$ f(x) = a \\sin(b x + c) $$\n",
    "\n",
    "Let's define the data and plot it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def target_function(x):\n",
    "    return # your function here\n",
    "\n",
    "def model_function(x, a, b, c):\n",
    "    return # your model here\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate some data\n",
    "num_points = 100\n",
    "noise_level = 1.0 # For modeling the experimental noise\n",
    "\n",
    "x = # your x values here\n",
    "data = # your data here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 203,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot the data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We indeed see some noisy oscillations, but we can still see the sinusoidal behavior. Our goal is to fit this data with the model function and extract the best-fit parameters.\n",
    "\n",
    "To do this, we need to minimize the residual between the data and the model function. The residual is usually defined as a difference between the data and the model prediction, scaled by the uncertainty in the data:\n",
    "\n",
    "$$ \\text{residual} = \\frac{y_i - f(x_i)}{\\sigma_i} $$\n",
    "\n",
    "where $y_i$ is the data, $f(x_i)$ is the model function, and $\\sigma_i$ is the uncertainty in the data.\n",
    "\n",
    "In the case of `lmfit`, it is usually convenient to provide all parameters of the model in a single dictionary. Let's define the model function and the residual function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def residual(params, x, data, uncertainty):\n",
    "    return # your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To create the initial guess of the parameters, we need to use the `lm.create_params()` function, and put the initial guess as arguments. Since we don't know anything the parameters yet, we can but $a=b=1$ for amplitude and frequency, and $c=0$ for the phase."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "params = # your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, once the params are created, we can use the `lm.minimize()` function to minimize the residuals w.r.t. the free parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "out = # your code here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The chi-square statistic is pretty large, which means that the fit is not optimal. Let's print the best-fit parameters and plot the data and the model function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot the data and the fit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We don't learn any thing... This is because the least-squares algorithm is sensitive to the initial guess of the parameters. Therefore, we need to provide better initial guesses. Let's try to do this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your code here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This time chi-square is much smaller, which means that the fit is better. Let's print the best-fit parameters and plot the data and the model function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot the data and the new fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This time our fit is perfectly fine. We can see that the best-fit parameters are very close to the true values."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "metatrain",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}